
# get the command arguments from perl
# rsm number
# runfolder
# run number

args=(commandArgs(TRUE))

for(i in 1:length(args)){
    eval(parse(text=args[[i]]))
}	 

rsm <- 4
runfolder <-"/home/rumen/simulation1/K=0_P=0_MU=1e-06_NU=1e-04_RSM=4_B=0.004_D=0.003_NN=100_T=14600_H=300_W=300" 
run <- 2

rsmdata <- read.csv("/home/rumen/cryptsim/analysis/run.simulation.coefficients.models.csv",header=FALSE)
fitness.log <- read.csv(paste(runfolder,"/rep",run,".fitness.log",sep=""),header=FALSE)

rsmmodel <- rsmdata[rsmdata[,1]==rsm,]
repr <- as.numeric(strsplit(as.character(rsmmodel[1,3])," ")[[1]])
surv <- strsplit(as.character(rsmmodel[2,3])," ")[[1]]
muta <- strsplit(as.character(rsmmodel[3,3])," ")[[1]]

repr.reduced <- repr[repr!=0]
relative.fitness.coefficients <- NULL 
for(i in 1:length(repr.reduced)){
 combos <- combn(repr.reduced, i)
 combos <- combos + 1
 relative.fitness.coefficients <- c(relative.fitness.coefficients, apply(combos, 2, prod))
} 
reproduction.fitness.classes <- c(1.00,sort(unique(relative.fitness.coefficients)))

# Remove the wild-type clone
fitness.log <- fitness.log[-which(fitness.log[,2]==0),]

# Prepare fitness landscape 4 by 8 = 32 cells
yy <- expand.grid(locus1 = c(0,1),  locus2 = c(0,1))
xx <- expand.grid(locus3 = c(0,1),  locus4 = c(0,1), locus5 = c(0,1))




combos <- digitsBase(0:31,base=2,n=5)
all.loci.combos <- (repr.reduced+1) * combos
all.loci.combos[which(all.loci.combos==0)] <- 1.0
fitness.combos <- apply(all.loci.combos,2,prod)
#orrd <- order(fitness.combos)
#combos <- combos[,orrd]
#all.loci.combos <- all.loci.combos[,orrd]
#fitness.combos <- fitness.combos[orrd]

# Get specific day
fitness.timepoint <- fitness.log[fitness.log[,1]==7300,]

genotypes.frequencies <- data.frame()
for(i in 1:nrow(fitness.timepoint)){
genotypes.frequencies[i,1] <- paste(fitness.timepoint[i,8:12], collapse="")
genotypes.frequencies[i,2] <- as.numeric(fitness.timepoint[i,3]) 
}
unique.genotypes.frequencies <- data.frame()
unique.genotypes.frequencies <- t(t(unique(genotypes.frequencies[,1])))
unique.genotypes.frequencies <- cbind(unique.genotypes.frequencies,rep(0,nrow(unique.genotypes.frequencies)))
for(i in 1:nrow(unique.genotypes.frequencies)){
      unique.genotypes.frequencies[i,2] <- sum(genotypes.frequencies[genotypes.frequencies[,1]==unique.genotypes.frequencies[i,1],2])
}

frequency.combo <- rep(0,ncol(combos))
for(i in 1:ncol(combos)){
      if(length(unique.genotypes.frequencies[unique.genotypes.frequencies[,1]==paste(combos[1:5,i], collapse=""),2])>0){
      frequency.combo[i] <- unique.genotypes.frequencies[unique.genotypes.frequencies[,1]==paste(combos[1:5,i], collapse=""),2]
      }
}
frequency.combo <- as.numeric(frequency.combo)

z <- matrix(fitness.combos,nrow=4,ncol=8,byrow=TRUE) 
colz <- matrix(frequency.combo,nrow=4,ncol=8,byrow=TRUE)


f1 <- c(0.3352588 ,-0.940788, -0.05019309   , 0)
f2 <- c( 0.4839025 , 0.126242 , 0.86596853  ,  0)
f3 <- c(-0.8083563, -0.314612 , 0.49757338  ,  0)
f4 <- c(0.0000000,  0.000000,  0.00000000  ,  1)
fov <- rbind(f1,f2,f3,f4)
open3d(userMatrix=fov,windowRect=c(0,0,350,250), viewport=c(0,0,350,250),zoom=1.1)
rg.ramp <- colorRampPalette(c("darkgreen", "yellow", "red","blue"))
# Make 10 shades of green to red, dont hurt reader's eyes
ncolors <- 10
col.arr <- rg.ramp(ncolors)
# Scale down by a factor so that the maximum height is 4, 0.02 is background division rate
wt.div.rate <- 0.02
newz <- (z*wt.div.rate) * (4.0/max(z*wt.div.rate))
zoffset <- wt.div.rate*(4.0/max(z*wt.div.rate))
shift <- 0.1
for(i in 1:nrow(z)){
for(j in 1:ncol(z)){
      # Make the 4 top points
      ttl <- c(i-1+shift,j-1+shift,newz[i,j])
      ttr <- c(i-1+shift,j-shift,newz[i,j])
      tbl <- c(i-shift,j-1+shift,newz[i,j])
      tbr <- c(i-shift,j-shift,newz[i,j])
      # Make the 4 bottom points 
      btl <- c(i-1+shift,j-1+shift,zoffset)
      btr <- c(i-1+shift,j-shift,zoffset)
      bbl <- c(i-shift,j-1+shift,zoffset)
      bbr <- c(i-shift,j-shift,zoffset)
      # Make 1 quad that is the top
      cover <- rbind(ttl,ttr,tbr,tbl)
      # Make the walls
      w1 <- rbind(btl,btr,ttr,ttl)
      w2 <- rbind(btr,bbr,tbr,ttr)
      w3 <- rbind(bbr,bbl,tbl,tbr)
      w4 <- rbind(bbl,btl,ttl,tbl)
      # Make walls transparent
      quads3d(rbind(w1,w2,w3,w4),col="lightgrey", line.antialias=TRUE)
      quads3d(cover, col=col.arr[ceiling((colz[i,j]+0.0000001)*length(col.arr))],lit=FALSE, line.antialias=TRUE,point.antialias=TRUE)
}
}
a2 <- format(2/4*max(z*wt.div.rate),digits=2)
a4 <- format(max(z*wt.div.rate),digits=2)
axis3d(c("z-+"),tick=FALSE, at=c(0,2,4),labels=c(wt.div.rate,a2,a4),line=0)
lines3d(rbind(c(0,8,2),c(4,8,2),c(4,0,2)),col="black",lit=FALSE)
# make bounding box comprised of 3 walls
# bottom plate
quads3d(rbind(c(0,0,0),c(4,0,0),c(4,8,0),c(0,8,0)),col="black", line.antialias=TRUE,front="lines")
# left plate
quads3d(rbind(c(0,8,0),c(4,8,0),c(4,8,4),c(0,8,4)),col="black", line.antialias=TRUE,front="lines")
# farther plate
quads3d(rbind(c(4,0,0),c(4,0,4),c(4,8,4),c(4,8,0)),col="black", line.antialias=TRUE,front="lines")
par3d(ignoreExtent=TRUE)
for(i in 1:ncolors){
      clr <- rbind(c(4-shift,-0.1,(i-1)*(4.0/ncolors)+zoffset),
      	           c(4-shift,-0.6,(i-1)*(4.0/ncolors)+zoffset),
		   c(4-shift,-0.6,(i-1)*(4.0/ncolors)+(4.0/ncolors)+zoffset), 
		   c(4-shift,-0.1,(i-1)*(4.0/ncolors)+(4.0/ncolors)+zoffset))
      quads3d(clr,col=col.arr[i],lit=FALSE, line.antialias=TRUE)
}
# 0, 50% 100 % frequeency
lines3d(rbind(c(4-shift,-0.6,ncolors*(4.0/ncolors)+zoffset),c(4-shift,-0.8,ncolors*(4.0/ncolors)+zoffset)),col="black",lit=FALSE)
lines3d(rbind(c(4-shift,-0.6,(ncolors/2)*(4.0/ncolors)+zoffset),c(4-shift,-0.8,(ncolors/2)*(4.0/ncolors)+zoffset)),col="black",lit=FALSE)
lines3d(rbind(c(4-shift,-0.6,zoffset),c(4-shift,-0.8,zoffset)),col="black",lit=FALSE)
text3d(c(4-shift,-1.5,zoffset), text="0%")
text3d(c(4-shift,-1.5,(ncolors/2)*(4.0/ncolors)+zoffset), text="50%")
text3d(c(4-shift,-1.5,ncolors*(4.0/ncolors)+zoffset), text="100%")

rgl.postscript("bla2.pdf","pdf")











# Add extra row and column
newz <- cbind(rep(0,4),z,rep(0,4))
newz <- rbind(rep(0,10),newz,rep(0,10))
newcolz <- cbind(rep(0,4),colz,rep(0,4))
newcolz <- rbind(rep(0,10),newcolz,rep(0,10))


#persp.withcol(x=c(1:4),y=c(1:8), heightz=z,z=colz,nb.col=30,theta=-60,phi=5)
persp.withcol(x=c(1:6),y=c(1:10), heightz=newz,z=newcolz,nb.col=30,theta=-60,phi=5)

persp.withcol <- function(x,y,heightz,z,nb.col,...,xlg=TRUE,ylg=TRUE) { 
colnames(z) <- y 
rownames(z) <- x 
nrz <- nrow(z) 
ncz <- ncol(z) 
color <- palette(rgb.palette(nb.col)) 
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz] 
facetcol <- cut(zfacet, nb.col) 
par(xlog=xlg,ylog=ylg) 
persp( as.numeric(rownames(z)), as.numeric(colnames(z)), as.matrix(heightz), col=color[facetcol], ... ) 
}



persp(z,col=colz,theta=-60,phi=5)


#points(2,2,pch=21,cex=1.5,col="red")

# For a given clone



# Get all the unique sampling time points
sampling.time <- unique(glog[,3])

# Set color palette
rgb.palette <- colorRampPalette(c("black", "white"))
palette(rgb.palette(max.fitness.category))




# PNG files to make a movie
dir.create(paste(runfolder,"/fitness.rep",run,sep=""))
for(i in 1:length(sampling.time)){
png(paste(runfolder,"/fitness.rep",run,"/rep",run,".",sprintf("%04d",i),".png",sep=""),width=900,height=900,units="px")
par(mar=c(0,0,0,0))
current.time <- sampling.time[i]
current.glog <- glog[glog[,3]==current.time,]
mat.selective.col <- matrix(final.fitness[current.glog[,5],3],nrow=300,ncol=300,byrow=TRUE)
color2D.matplot(mat,cellcolors=mat.selective.col,xlab=NA,ylab=NA,do.hex=TRUE,border=mat.selective.col,axes=FALSE)
dev.off()
}

# mencoder mf://*.png -mf w=900:h=900:fps=2:type=png -ovc copy -oac copy -o output.avi
